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ABSTRACT 

In this paper the classical Euler-Bernoulli beam (CEBB) theory is reformulated utilising fractional 
calculus. Such generalisation is called fractional Euler-Bernoulli beams (FEBB) and results in non-local 
spatial description. The parameters of the model are identified based on AFM experiments concerning 
bending rigidities of micro-beams made of the polymer SU-8. In experiments both force as well as 
deflection data were recorded revealing significant size effect with respect to outer dimensions of the 
specimens. Special attention is also focused on the proper numerical solution of obtained fractional 
differential equation. 


1 Introduction 

A quantitative understanding of a size effect in engineering materials is of great importance in the 
design phase of micro- and nano-electromechanical systems (MEMS and NEMS), and as a technique 
of homogenisation regarding materials with real micro-structure. A size effect is reflected, for example, 
in a stiffer elastic response to external forces on small scales. This has been recognised in several 
experiments on metals and polymers, for example on copper (Fleck et al. 1994 [17]), silver (Ma 
and Clarke 1995 [24]), zinc oxide (Stan et al. 2007 [35]), lead (Cuenot et al. 2004 [12]), carbon 
nanotubes (Salvetat et al. 1999 [33]), epoxy (Chong 2002 [10]) and polypropylene (McFarland 
and COLTON 2005 [25]). Actually, a physical reasoning for the origin of a size-dependent material 
behaviour, such as non-negligible complex interactions of molecular chains in polymers (NlKOLOV et 
al. 2007 [27]), the rearrangement of atoms or molecules near the surface, the influence of the grain-size 
in polycystals (Smyshlyaev and Fleck 1996 [34]) and long-ranging influences from dislocations, 
voids or some inner micro- or nano-structure, is manifold and poorly understood especially in combined 
effects. For the reason that conventional continuum theory is unable to predict size effect, different 
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non-conventional continua are proposed in literature, like non-local theories (PEDDIESON et al. 2003 
[31], Eringen 2010 [16]), strain-gradient theories (TOUPIN 1962 [39], Mindlin and Eshel 1968 
[26]), micropolar theories (Eringen 1966 [15], Nowacki 1972 [28]) or theories of material surfaces 
(Gurtin & Murdoch 1975 [18]). 

In addition to these, the fractional calculus is assumed to be a promising candidate for modelling 
scale-dependent material behaviour too, by using a concept of intrinsic fractal structures to describe a 
potential disordered material’s micro-structure (cf. Carpinteri 1994 [8]). For detailed descriptions 
see Klimek (2001) [20], Vazaquez (2004) [41], Lazopoulos (2006) [21], di Paola et al. (2009) 
[13], Atanackovic & Stankovic (2009) [1], Carpinteri et el. (2011) [9] or Drapaca & Sival- 
OGANATHAN (2012) [14]. The formulation of fractional elasticity stated here defines spatial derivatives 
of arbitrary order including the positive features - in contrast to the specified non-conventional con¬ 
tinua - (1) to introduce a smaller number of additional material parameters, (2) to fit into the general 
framework of classical continuum mechanics and (3) to possess clear physical interpretation (e.g. re¬ 
garding to the intrinsic length scale). 

In Sect. 2.1 , we start with a description of a generalised continuum mechanical framework including the 
specific formalism of fractional calculus proposed in Sumelka 2014 [36]. More precisely, we define 
the gradient of the motion of a body, i.e. the deformation gradient, with the help of an differentiation 
operator a of real order. Having this, strain and stress measures, as well as a formulation of the con¬ 
servation of momentum will be derived according to conventional kinematics of small deformations. In 
Sect. 2.2 the classical Euler-Bernoulli assumptions will be utilised to define a displacement field 
valid for bending of slender beams. By applying equilibrium considerations to an infinitesimal beam 
element, the fractional Euler-Bernoulli differential equation will be presented as the basis for the 
numerical study. Section 3.1 clarifies the boundary conditions for the model in accordance with the 
experimental realisation. We will investigate the bending behaviour of micro-beams with rectangular 
cross-sections that are clamped on one side and loaded at the other by a single concentrated force. 
The computational algorithm, explained in more detail in Sect. 3.2, will use a one-dimensional mesh 
of equidistantly distributed calculation points along the beam’s axis. In between these nodes, the frac¬ 
tional Euler-Bernoulli differential equation will be solved including additional conditions for the 
transition of moments and forces between the sub-elements of the model. However, the operator for the 
differentiation of real order requires an additional step of extended integration over the neighbourhood 
of the relevant point to collect and weigh the non-local information of the body. The step of integration 
will be implemented using the fractional trapezoidal scheme Odibat (2006) [29], Blaszczyk et al. 
(2013) [6] and SUMELKA & BLASZCZYK (2014) [37], All these preparations allow us to build a stiff¬ 
ness matrix of the system, to be solved by the LU decomposition method. The benchmark numerical 
results will be presented in Sect. 3.3 for beams of different (thicknesses/lengths) and different fractional 
parameters will be applied to optimise the results to fit to the experiments. 

Finally, in Sect. 4 with the help of an Atomic Force Microscope (AFM) force as well as deflection data 
of micro-beams made of the polymer material SU-8 the fractional beam model will be identified. The 
experimental approach and the measured data will be explained to emphasise that the bending stiffness 
of the SU-8 structures revealed a size effect on the thickness and the length of the samples. 
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2 Fractional Euler-Bernoulli beams 


2.1 Small strain fractional elasticity 

There are many ways to generalise continuum mechanics in terms of fractional calculus, as mentioned 
in Klimek (2001) [20], Vazaquez (2004) [41], Lazopoulos (2006) [21], di Paola et al. (2009) 
[13], Atanackovic & Stankovic (2009) [1], Carpinteri et el. (2011) [9] or Drapaca & Sival- 
OGANATHAN (2012) [14]. Herein we follow the concept presented in Sumelka (2014) [36]. Because 
the FEBB theory results from putting specific restrictions on general small strain fractional elasticity, 
we start with a short introduction to the latter one - for details cf. Sumelka (2014) [36], Sumelka & 
Blaszczyk (2014) [37]. 

The concept of fractional continua discussed in Sumelka (2014) [36] is based on fractional deforma¬ 
tion gradients concept, namely 

F(X,i) =£«~ 1 D a ( j){X,t), or F aA = t^ 1 D °0 a e a 0 E A , (1) 

XX X Xa 

and 

F(x, t) = l^T x D a p>(x, t), or F Aa = £%~ 1 D a ip A F, A <8> e a , (2) 

XX X Xa 

where F and F arc fractional deformation gradients in material and spatial descriptions, respectively, 

x x 

lx and l x are corresponding length scales, cf) defines the regular motion of the material body while 
ip its inverse. The fractional differential operator D a is defined as the Riesz-Caputo (RC) fractional 
derivative 

R lA/(t) = (JWM + (- 1 )" tO?m) , O) 

where a > 0 denotes the real order of the derivative, D denotes ’derivative’ (RC stands for RlESZ- 
Caputo), a,t,b are so called terminals, T is the Euler gamma function, ^Dff{t), are 

left and right Caputo’s fractional derivatives, respectively PODLUBNY (1999) [32], KlLBAS ET AL. 
(2006) [19], LESZCZYNSKI (2011) [22]. From the definition of the RC derivative by Eq. (3) it is clear, 

that fractional deformation gradients given by Eqs (1) and (2) are non-local (each component of F and 

X 

F governs the information from the surrounding, described by terminals a and b). In consequence, other 

X 

related measures of deformation (e.g. strain tensor) which are obtained through the relations analogous 
as in classical case, are also non-local. To conclude, fractional kinematics is controlled by new material 
parameters: length scale, order of derivation and type of fractional derivative Oliveira & Machado 
(2014) [30] - they should be identified for specific material. 

Thus, the finite fractional strains are obtained from the difference in scalar products in actual and refer¬ 
ence configurations (cf. Fig. 1) 


E = -(F 
2 k 


O ^ 

t F-I), 


1,0 0 

or Eab = ^{F Aa F aB — Iab)^A®~Eb, 


1 O T <> 1 o o 

e = 2^” F ~ F_1 )> ° r eab = ^ab- F aA F Ab) e a® e b i 


(4) 

(5) 


where E is the classical Green-Lagrange strain tensor or its fractional counterpart (symmetric), e 

o 

is the classical Euler-Almansi strain tensor or its fractional counterpart (symmetric), and F can be 

_ OL 

replaced with F or F or F or F (cf. Fig. 1). 


3 




Wojciech Sumelka, Tomasz Blaszczyk, Christian Liebold - Compiled on 2015/02/06 at 01:34:59 



Figure 1: The relationship between material and spatial line elements with their fractional counterparts 
F = FF _1 F F = FF and F = FF _1 . 

X x x x XX 


An infinitesimal fractional strain is obtained as in the classical set-up by considering the relationship 
between fractional strains and fractional displacement gradients tensors and omitting higher order terms 
in obtained relations. Hence, infinitesimal fractional Cauchy strain tensor is defined as (assuming that 

^ = f-x = 4 ) 


GradU + GradU T 

1 

gradu + gradu 7 

X X 

2 

X X . 


( 6 ) 


where 

GradU = F - I, 

or r~ 1 D a U a = 

- (-^ aA laAij^a ® ^^4 


X X 

X A 

X 

and 

gradu = i — F, 

or r~ 1 D a u A = 

{i'Aa F Aa)^A ® 


XX Xa x 


In the above equations U and u denote material and spatial displacement. 

It is clear that for the fractional kinematics appropriate fractional kinetics should be defined SUMELKA 
et al. (2014) [38]. It can be shown that fractional stress field a satisfies analogous relations as in the 
classical formulation. Assuming that the conservation of mass holds, hence 

Pq&V = fJ()dV = pdv = pdv, (9) 


or shortly 

o 

Po = Jp, (10) 

where po (po) is the reference mass density (fractional counterpart), p (p) is spatial mass density (frac- 

o o 

tional counterpart), and J = detF is a jACOBian, we have 
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where v is a velocity, and f is a body force per unit mass, so the local form is 



or in the absence of inertia forces 

div<x T + pf = 0. (13) 

Furthermore, the symmetry of a (so a or a) can be checked in a classical manner utilising the balance 
of moment of momentum. The result given by Eq. (12) gives justification and is analogical for the one 
postulated in Atanackovic & STANKOVIC (2009) [1] for similar small strain elasticity generalisation. 

Finally, the small strain fractional elasticity is governed by (for purpose the denotation (•) is omitted) 

&ij,j T bj = 0. 

(D a Ui + D a u-j 

J 2 \Xj Xi J 

< O'v = c ijki £ kh (14) 

Ui = Ui , X E Slu, 

<JijTlj — ti, X £ Her; 
i fl( 7 nn cr = 0 and fij/U fio-= fi. 

In the above equation we have denoted: b is the body force, C e is the stiffness tensor, f Ijj and Q a are 
parts of boundary 11 where the displacements and the tractions are applied, respectively. It is clear that 
the classical (local) solution appears as a special case of the proposed generalisation when a = 1. 

2.2 Governing equations 

Governing equations for the FEBB are obtained under the analogous assumptions as in a classical case 
(for classical case cf. Bauchau & CRAIG (2009) [3]). Without loss of generality we consider the effect 
of extension and bending (other concepts like Saint-Venant assumptions associated with torsion can 
be incorporated by analogy). 

The classical ad hoc assumptions for kinematics hold, namely: 

Assumption 1 The cross-section is infinitely rigid in its own plane, 

Assumption 2 The cross-section of a beam remains plane after deformation, 

Assumption 3 The cross-section remains normal to the deformed axis of the beam. 

Within the above assumptions we can express the 3D displacements in terms of the ID beam displace¬ 
ments (direction 1 coincides with a beam axis): 

ui(x 1 ,x 2 ,x 3 ) = tti(xi) + tc 3 $ 2 (a:i) - 

u 2 (x 1 ,x 2 ,x 3 ) = u 2 (x i), (15) 

u 3 (x 1 ,x 2 ,x 3 ) = u 3 (x i). 

In Eqs (15) we have assumed the following sign convention: the rigid body translations of the cross- 
section U) (x'i), v, 2 {x \), and /^(aq) are positive in the positive direction of the coordinate axis, respec¬ 
tively; the rigid body rotations of the cross-section ^(aq) = — d ^^ 1 ' 1 and ^(aq) = are 

positive keeping the right-hand rule. 


_ Wa-l 
— 2 *" 


u ij + u j t i 
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Because finally the experimental set-up considers bending in (l,3)-plane we reduce the set of Eqs (15) 
to 

dtt 3 (xi) 

Ul(X 1 ,X2,X 3 ) = -X 3 — -, 

dxi 

U2{x l,X2,X 3 ) = 0, (16) 

U 3 (xi,X 2 ,X 3 ) = U 3 (xi). 

Applying the displacement relations Eqs (16) to the definition of fractional Cauchy strain Eq. (6) we 
observe that the only non-zero element is 

til = ~x 3 t^ x D a { ^ U ^ Xl \ = zl a ~ l k 2 = -zi a ~ x w, x i , (17) 

x, dx\ 

so 

(Til = -zt^Ew^i , (18) 

where z = x 3 , E is Young’s modulus, w = u 3 (xi) and « 2 (xi) = — w,-^. It is clear that the transverse 
shear stress due to flexure cannot be obtained in this way due to the third Euler-Bernoulli assump¬ 
tion - it is determined from equilibrium considerations as in the classical set-up (we further denote it by 

< 3 )- 



Figure 2: The sectional stress resultants. 


Next, following Fig. 2 the kinetic variables, called the sectional stress resultants, are 

M 2 = [ (TiiXsdA = — w n i f z 2 ^- 1 EdA, (19) 

Ja J a 

and 

F 3 = f (7* 3 dA. (20) 

JA 

Up to this point we have three unknowns: displacement u 3 , (fractional) strain «2 of the beam, and stress 
resultant M\. Thus far, we have obtained two equations which combine them, namely: 

• strain-displacement relation (12(^1) = —-P a ( d ^^ 1 ) ), and 

• constitutive relation M 2 = K 2 J A z 2 £ a ~ 1 EdA. 
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The missing last one is derived using equilibrium considerations, so called NEWTONian method. 



Figure 3: Free body diagram for the transverse shear forces and bending moments. 


Thus, in the absence of distributed moments through the beam axis, following the free body diagram 
shown in Fig. (3), we have the following relations 


and 


=-»<**>• 

( 21 ) 

0 . 

ax\ 

( 22 ) 


where ps(xi) denotes the distributed load. Next, by taking a derivative of Eq. (22) then introducing 
Eq. (21), the bending moment equilibrium equation is 


d 2 ^2 , x 

—^2“ = -P3{Xl), (23) 

which defines the third lacking equation, for the full description of beam bending in the plane (1,3) (in 
the absence of distributed moments through the beam axis). 

Finally, combining Eq. (19) and Eq. (23) we obtain a single equation in terms of w 


e- ] Ehu niu = 


where / = J A z 2 dA. When order of fractional continua 


P3(xi), (24) 

is taken a = 1 classical solution is recovered 


£°EIw, 1 in = El 


d 4 w 

dxf 


P3(xi)- 


(25) 


3 Numerical study 

3.1 Bending test 

The real physical conditions in micro-beam bending tests, discussed in Sect. 4, correspond to a can¬ 
tilever beam model, loaded at the free end by a point load P - cf. Fig. 4. 
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Figure 4: Computational model of a micro-beam bending test. 


For such a configuration resultant moment is known, namely 

M 2 (x t) = PL - Px i, 


thus Eq. (23) is satisfied. Therefore the problem is governed by (cf. Eq. (19)) 

'-D a {^-)r~ 1 EI = PL — Px 1 , 

Al 1 

< w(x i = 0) = 0, 

.&(*!= 0 ) = 0 , 

where / = • 


(26) 


(27) 


3.2 Computational algorithm 

In this section we present a numerical scheme for the problem (27). We introduce the homogeneous 
grid of nodes (see Fig. 5). We see that additional fictitious nodes x]~ m ,..., aq 1 and x ^ +1 ,..., x^ +m 
placed outside the domain [x®, x±] are introduced. We denote a value of the beam deflection at the node 
x\ as w (x \) = w t . Next, by analogy as in ClESlELSKl & LESZCZYNSKI (2006) [11], Sumelka & 
Blaszczyk (2014) [37], we assume that for all fictitious nodes on the left the beam deflections are 
w- m = w - m+ 1 = ... = w -1 = w o and on the right the beam deflections are wj\r+ m = ww+m-i = 
... = w N+ 1 = w N . 


L =NAx , 


x 


-m 

1 








< - 2 1 -► 




N -1 
1 


X 


N 

1 


X 


N + m 
1 


Figure 5: Spatial discretization for micro-beam. 

The problem of numerical solutions of equations containing simultaneously the left and right frac¬ 
tional derivatives has recently been drawing the attention of many authors Asad et al. (2014) [2], 
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Blaszczyk & Ciesielski (2014) [5], Blaszczyk (2015) [4], CRESSON et al. (2013) [7], Xu & 
Agrawal (2014) [42]. In this paper we use a scheme which is based on the fractional trapezoidal rule 
Odibat (2006) [29], Blaszczyk et al. (2013) [6] and Sumelka & Blaszczyk (2014) [37]. The 
left Caputo derivatives at nodes x\ is approximated by the formula 



dm 
dxi ^ 


x\=x\ 


1 f w" (r) dr 

r(l - a) J (x\ - r)“ 


r i—m 
V 1 


(An) 


1—a 


T (3 — a) 

i— 1 


| (m — l) 2 a — (m + a — 2) m 1 “ w" (x\ m ) + w" (x\) 


+ E [(* - J - i) 2 ' 0 - 2 (* " + (i - J + 1) 

j=i—m +1 


2—a 


w" X J 


0} 


(28) 


where (•)" denotes a second derivative with respect to x\. Using the central difference formula for 
the second order derivatives occurring in (28) we finally obtain the approximation of the left Caputo 
derivative 


'. D a w 

i — m 1 ^X\ u_/ 


(Axi) 1 
T (3 — a) 


i +1 

E w 3 v l( i ’j)i 

j—i—m— 1 


(29) 


where 


v\ (i,j) 


(m — l) 2 “ — (m + a — 2) m 1 a 

for j = i — to — 1 

(m - 2 f~ a - 4 (m - 1 f~ a 
+2 (m + a — 2) m l ~ a + m 2 ~ a 

for j = i — to 

-2 (m - 2) 2- " + 5(m - 1) 2_ “ 

— (m + a — 2) m l ~ a — 2 m 2 ~ a + 1 

for j = i — to+1Ato = 2 

(to - 3) 2- " - 4(to - 2) 2_q + 6(m - 1) 2 "“ 

— (m + a — 2) to 1- ® — ?n 2 ”" 

for j = i — m+lAm>2 

(* - j ~ 2) 2_a - 4(i - j - l) 2 “ a + 6 (i - jf~ a 
-4(i-j + l) 2 ~ a + (i-j + 2) 2 - a 

for j = i — to +2, ! 

6 - 4 • 2 2 " q - 3 2-q 

for j = i — lAm>2 

2 2_ “ — 4 

for j = i 

1 

for j = i + 1 

0 

otherwise 
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We determine the discrete form of the right fractional CAPUTO derivative in a similar way 


c no / d- w \ 
11 


i+m 
L 1 


— 1 f w" (r) dr 

r(l-a) J (t — x\) a 


- (Axi) 


1 —a 


r (3 — a) 

z+m+1 


(m — l) 2 “ — (m + a — 2) m 1 a w" (+" m ) + w" (a+ 


+ Y ti-i~ 1) 2 "“ - 2 (j - if~ a + (j-i + l) 2 -“ 


j=i-l 

- (Ax,)" 1 " 0 i+ ^t‘ 
r (3 - a) A (.,]), 

j=l~ 1 


W X 


(31) 


where 


V2 (i,j) = < 


(m — l) 2 “ — (m + a — 2) m 1 a 

for j 

= i + m + 1 

(m - 2) 2 "“ - 4 (m - 1) 2_Q 



+2 (m + a — 2) m 1- " + m 2 ~ a 

for j 

= i + m 

-2 (m - 2) 2 ~ a + 5(m - 1) 2_ “ 



— (m + a — 2) m 1_a — 2 m 2 ~ a + 1 

for j 

= i + ra-lAm = 2 

(m - 3 ) 2 " q - 4(m - 2) 2 " Q + 6(m - 1) 2 “" 



— (m + a — 2) m 1-a — ?n 2_ “ 

for j 

= i + ra-lAm>2 

(j-i- 2?~ a - 4(j - i - l) 2 ~ a + 6 (j - if~ a 



— 4(j — * + l) 2_a + (j — i + 2) 2 ~ a 

for j 

= i + 2,..., i + m — 

6 - 4 • 2 2-a - 3 2_a 

for j 

= j + lAm>2 

2 2_ “ — 4 

for j 

= i 

1 

for j 

= i — 1 

0 

otherwise 


(32) 

Finally, we present the discrete form of the considered problem (27). For calculation of values 
wq, w\,..., wn we need to solve the system of A r + 1 linear equations. For every grid node x\, where 
/ = 0,..., N, we can write the following equations 


tuo = 0, 

wi = 0 , 

< z+1 z+m+1 

< j=i—m —1 j=i— 1 


P (L — zAxi) 


(33) 


where 

l a ~ x EI 

(Axi) 1+Q r (3 — a) ' 


(34) 


The parameter m, in above formulae, represents the number of spatial nodes for which the fractional 
derivatives, Eq. (28) and (31), are calculated. The relation between the parameter m, length scale i and 
step size Aaq is expressed by: Aaq = + Thus, for larger values of t, rri should also be increased 
appropriately, because the convergence of solution, represented in Eq. (33), is in relation to Aaq. 
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3.3 Benchmark results 

As in the CEBB theory, also in its fractional non-local version, the information about the cross section 
is introduced through the moment of inertia I only. In this sense, the FEBB theory makes it possible to 
analyse the non-local effects through beam length L. In the following series of benchmark results we 
consider the influence of length scale £ and order of fractional continua a on cantilever beam bending 
test results. 

On the basis of the numerical scheme presented in the previous section we implemented a computer 
program in Maple and carried out computational simulations for different values of parameters a, £ and 
m. We applied the LU decomposition method in order to numerically solve the system of equations (33). 
In all presented benchmark examples we assumed P = 1 N, E = 1 Pa, L = lm, W = T= ^L 
and Axi = 0.001 m. We denote by w (L, a) a value of a beam deflection at the loaded end w (xf ) for 
fixed a £ {0.4, 0.5,0.6, 0.7, 0.8,0.9,0.99} and by w (L, 1) a value of a beam deflection at the loaded 
end w (xi) for classical local i.e. a = 1. In Fig. 6 the numerical results are presented. Two cases are 
considered: (i) the interval of fractional differentiation 21 is smaller (or equal) than the smallest beam 
dimension T, and (ii) 21 is smaller (or equal) then beam length L. 



Figure 6: The relation between the normalised cantilever beam deflection at loaded end vs. order of 
fractional continua a for different length scales /(Left: case (i). Right: case (ii). 


We observe that when length scale £ increases in comparison to beam length L, the non-local effects are 
more vivid. If £ approaches zero and/or a one, the FEBB reduces to the CEBB model. It is also clearly 
seen that for smaller values of a the difference between a classical and fractional result increases. 
We can conclude that both £ and a control the stiffness of the fractional beam - this result is crucial 
concerning model identification in the next Sect. 4. 


4 Experimental validation 

4.1 Experimental set-up 

Deflection and force data from real micro-beam bending tests were recorded with an off-axis laser- 
reflective Multiview-1000 AFM-stage from Nanonics Imaging Etd. The system is composed of a flat 
scanner, consisting of a fine thread driven by piezo-elements with high-voltage power supply and a 
detection device working with four Photo-Sensitive Diods (PSD) interconnected as a WHEATSTONE 
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bridge to monitor deflections of a laser-beam path. The laser is reflected in an obtuse angle from a fixed 
AFM-cantilever such that the system directly monitors its deflections w c , when deformed by the piezo 
uplift (referred to as separation z). With the knowledge of the well calibrated spring constant of the 
AFM-cantilever of k c = 31.4 N m 1 , we were able to convert the PSD-signal into force P = k c w c (cf., 
Fig. 7). The calibration process used here is described in more detail in Varenberg et al. (2005) [40] 
and was performed using a precise silicon normal that was provided by the PTB (Physikalisch Tech- 
nische Bundesanstalt - Braunschweig). The AFM-cantilever is a custom-built cylindrical and curved 
cantilever made of glass, having a tip radius of about 20 nm. In micro-beam bending tests it can be as¬ 
sumed that the pure AFM-data consist of a combined signal of the deflection of the AFM-cantilever and 
the micro-beam’s deflection w in the following manner: z = w + w c . Hence, we were able to separate 
the deflection data of a micro-beam from the pure AFM-data in combination with the corresponding 
load, cf. Fig. 7. 



—measurement on 
rigid surface (I) 
- -measurement on 
micro-beam (II) 


80 
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,50 
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-600 -400 


-200 0 200 400 600 800 1000 1200 1400 

separation [nm] 


-600 -400 -200 0 200 400 600 800 1000 1200 1400 

separation, w [nm] 


Figure 7: Some exemplary original PSD-data (left) as well as shifted and converted data (right). 


The ratio P/w is called bending stiffness with the classical relationship to the elastic modulus E, width 
W , length L and thickness T 


4L 3 P 
WT :i w ’ 


(35) 


when the CEBB theory and a rectangular cross-section of the beam are assumed. The length and width 
of the samples have been measured in an optical microscope with a magnification of 500 times, whereas 
the thickness was taken to be the mean value from two different optical determination systems and 
values of a Scanning Electron Microscope (SEM). Thereby, the lengths were specified between the 
fixation of the cantilever on the solid glass support and the force application point of the AFM-tip, 
which could be varied. 

SU-8 is known as the photo-resist Nano™-SU-8 of the company MicroChem, used in the micro-system 
technology as described in Lorenz et al. (1997) [23]. Manufacturing of the samples was earned out in 
the following steps: 


• A 4-inch silicon wafer was coated with a thin metal film as a barrier layer. 


12 










Wojciech Sumelka, Tomasz Blaszczyk, Christian Liebold - Compiled on 2015/02/06 at 01:34:59 



Figure 8: An exemplary SEM picture of the free end of a tested micro-beam made of SU-8. 


• First, the viscous SU-8 resist was dissolved in solvents. A conventional spin-coating machine 
was applied to disperse the resist on the metal-coated silicon wafer. The rotational speed of the 
coating machine determined the thickness of the homogeneous layer (between 8—40 microns). 

• The solvent evaporated in the rotary process in large part and the remaining part was evaporated 
in a subsequent drying process at temperatures of 60°C - 95°C, whereby the material finally 
received its rigidity. 

• Structuring was achieved by a Laser Direct Imager (LDI). After exposure to light an additional 
heat treatment was earned out at 60°C - 95°C to assist the chemical reaction of illumination. 

• By using a proper developer, the regions of exposed SU-8 were dissolved from unexposed re¬ 
gions. 

• Followed by a chemical etch process, which does not attack the SU-8, the thin metal film on the 
silicon wafer was dissolved and the micro-beam structures were finally peeled off. 

• In a last step the micro-beam structures were glued on a support made of glass and fixed on an 
AFM sample holder. 

Samples of different lengths and widths have been produced in different thicknesses (T = 8.4 pm and 
14.4 pm) at the Fraunhofer Institute for Reliability and Microintegration, Berlin, and were tested with 
the AFM-technique according to Eq. 35. 
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Figure 9: The effect of the length of a beam made of SU-8 on the elastic modulus recorded with the 
AFM-technique. 


Results show that the values of elastic moduli E increase when the thickness of the cantilevers is 
decreased. The increase of elastic modulus correlates positively with the length L (see Fig. 9) and 
the values have nearly doubled in the case of the thinnest cantilevers. The size effect, that has been 
detected here, manifests itself in the slopes of the partial regression lines of the experimental data 
which show an overall increase with decreasing beam thicknesses. On the basis of general predictions 
of the material surface theory (Gurtin & MURDOCH 1975 [18]), regression lines were constructed as 
inversely proportional functions and independent of the actual thickness. We observe a dependency of 
E on the length L as well as on the thickness T. 

4.2 Model identification 

The discussed experimental results show a very complex behaviour of the real micro-beam. The scale 
effect is vivid when dimensions of rectangular prism specimen (L , W. T ) are changed, and is more 
pronounced when these dimensions arc comparable with an assumed real intrinsic length scale. In this 
sense, from the point of view of modelling, a 2D or even 3D description should be applied. Additionally, 
anisotropic length scale effects could be covered in a 2D or 3D description. Nevertheless, we claim that 
the FEBB model described here is able to model the present size effects with respect to a change of beam 
lengths and beam thicknesses by finding appropriate combinations of the two additional introduced 
parameters a and l, together with the type of fractional derivative (in our case the Riesz-Caputo 
one). However, similar to the CEBB model, the information about the cross section is only included 
in the second moment of inertia (cf. Eq. 27). Thus, from the first sight one can expect, that only scale 
effects between different micro-beam's lengths (L) can be covered by FEBB, but not between different 
thicknesses or widths (T or W). Fortunately, such a statement is not true at all. In our observations, a 
change of both, length and cross sectional area of micro-beams can be covered by FEBB formulation, 
by means of changing the order of fractional continuum 0 < a < 1. In Fig. 10 the comparison of the 
experimental results (cf. Fig. 9) and the results of the FEBB model are presented. Non-local moduli 
Enl have been calculated according to Eq. 35, utilising numerical deflections w (L, a). We observe, 
that the length scale in our FEBB model should be i = 60 pm, whereas Young’s modulus should be 
E = 6.9 GPa. The effect of a change of cross sections is modelled using a different order of fractional 
continua, namely: for T = 8.4 pm and W = 82 pm we have a = 0.8; and for T = 14.4 pm and 
W = 122 pm we have a = 0.4. It should also be emphasised, that the experimental result, that when 
the length L of micro-beams is growing, the non-local elastic modulus E\,\ converges to a single value. 
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This gives evidence to the previously expected behaviour of the FEBB model, that if the dimensions 
of beams become significantly larger than the length scale ( L , I f 7 . T A> £), the non-local modulus E'nl 
converges to the classical elastic modulus E. 



5 Conclusions 

In the present work, the non-local fractional Euler-Bemoulli beam theory is formulated as a generali¬ 
sation of classical Euler-Bernoulli beams, utilising fractional calculus. The model is implemented 
using codes developed by the authors and identified based on AFM experiments concerning bending 
of micro-beams made of the polymer SU-8, where scale effect is severely manifested. The fractional 
Euler-Bernoulli beam model introduces two additional parameters a and i together with the type 
of fractional derivative (in the discussed case the Riesz-Caputo one). It is shown, that this new model 
is able to give, in a qualitative and quantitative manner, good approximation of the revealed experimen¬ 
tal results. 
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